Stress Response in Confined Arrays of Frictional and Frictionless Particles 
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Stress transmission inside three dimensional granular packings is investigated using computer 
simulations. Localized force perturbation techniques are implemented for frictionless and frictional 
shallow, ordered, granular arrays confined by solid boundaries for a range of system sizes. Stress 
response profiles for frictional packings agree well with the predictions for the semi-infinite half 
plane of classical isotropic elasticity theory down to boxes of linear dimensions of approximately 
forty particle diameters and over several orders of magnitude in the applied force. The response 
profiles for frictionless packings exhibit a transitional regime to strongly anisotropic features with 
increasing box size. The differences between the nature of the stress response are shown to be 
characterized by very different particle displacement fields. 

PACS numbers: 62.20.-x 81.70.Bt 83.80.Fg 
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I. INTRODUCTION 

The manner in which granular materials respond to 
external loads belongs to a class of problems related to 
mechanical stability. Questions surrounding mechanical 
stability of particulate media have broad economical, en- 
vironmental, and societal applications. These include the 
common practical problem of a collapsing grain silo, nu- 
merous natural events such as avalanches and mud slides, 
along with other unexpected destructive occurrences like 
dam breakage. Many such events are reported annually 
and readily convey the dramatic, and often catastrophic, 
effects that are primarily due to the failure of various 
granular systems. 

The general physical theme that connects such a dis- 
parate range of problems is that of determining the un- 
derlying stress state of the system. Given information 
on the distribution of stresses within the system, can we 
predict whether such systems will remain mechanically 
stable, and what are the conditions that might promote 
or inhibit failure? One possibility is to employ estab- 
lished theories such as continuum, linear, elasticity the- 
ory Although elasticity theory applies over a wide 
range of situations, deviations from classical elasticity 
theory show up on a regular basis as some of the assump- 
tions that underpin the theory break down. These prob- 
lems have largely been emphasized over the past 15 years 
or so in studies of packings of granular materials. The 
essential point to note is that granularity, or the discrete- 
ness of the particles that constitute the packing, play a 
vital role in the resulting stress characteristics of the sys- 
tem due to a lack of a separation of length scales. Thus, 
it is becoming increasingly necessary to probe the rele- 
vant limits of classical elasticity theory and its range of 
valid application. For instance, classical elasticity is often 
assumed to remain valid not only at macroscopic scales 
to infer stress properties from displacement fields during 
loading events of the lithosphere [2], but also at the mi- 
croscopic level of interrogation in studies of nanoscopic 
indentation in relation to physical metallurgy [sl, 0] • 

A convenient method to test the range of validity of 



classical elasticity theory in granular systems that has 
been utilized over the past decade is the localized force 
perturbation technique. This procedure has the particu- 
lar benefit that for an isotropic elastic material the stress 
profile in response to localized forcing - Green function 
response - can be determined without the need for any 
free parameters. We review this result for the case of 
three dimensions - Boussinesq equations - in the follow- 
ing section (and see the Appendix). To summarize efforts 
to date, studies on two- and three-dimensional granular 
packings indicate that disordered and frictional systems 
tend to follow the predictions of isotropic elasticity the- 
ory. Whereas, ordered packings composed of frictionless 
particles result in strongly anisotropic stress profiles that 
are not consistent with a linear, isotropic, elastic model. 
Such anisotropic profiles (see below for a more qualitative 
definition) have been debated in relation to anisotropic 
elasticity theory based on differential equations of the el- 
liptic class ^ , and alternative descriptions of the hyper- 
bolic wave variety [l, • To distinguish between the tra- 
ditional isotropic result and anisotropic stress profiles we 
refer to these as one-peak and two-peak response profiles 
respectively, in reference to two dimensional experiments 
and simulations that have reported these observations 

In recent computer simulations the role of structure 
for frictionless granular packings on the stress response 
in a pseudo-infinite medium was investigated ^ . These 
studies showed that for ordered and quasi-ordered ar- 
rangements of frictionless grains, the response function 
was two-peak but consistent with the framework of an 
anisotropic elastic theory. A phenomenological model 
was employed to fit these anisotropic profiles that nec- 
essarily captured the two-peak character of the response 
profiles. As the structural disorder of the packing was 
increased, the stress response crossed over to a one-peak 
response, more in line with isotropic behavior. 

The problem addressed here is to determine how con- 
finement affects the stress response of granular packings. 
In most previous studies samples were used of sufficient 
size that the perturbed packings could be considered as 
semi-infinite in extent. In some cases the roughness of the 
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supporting base in packings prepared under gravity was 
seen to influence the magnitude of the stress response, 
although the response profiles retained isotropic, elastic 
character [l6|. We also highlight two particularly rele- 
vant efforts that have been proposed to describe stress 
response specific to two dimensional systems. Exten- 
sions to three dimensions and especially for fully con- 
fined, shallow packings remains an ongoing effort. A 
non-linear elasticity formalism was developed [I^ that 
captures the main features of experimental results for 
two dimensional granular packings [8[. The isotropic re- 
sult is recovered through a multipole expansion of the 
normal stress response component cFzz of classical elas- 
ticity. This theory was also able to describe anisotropic 
profiles by including generalized nonlinear strain com- 
ponents through the introduction of phenomenologically 
defined elastic constants. A future extension of this par- 
ticular formalism may be useful in describing the stress 
response for a wide range of granular systems in both 
two and three dimensions. Additionally, force perturba- 
tion simulations of confined Lennard- Jones glasses [HI 
found a range of stress response profiles that appear, at 
least superficially, similar to some of the isotropic results 
we present here. Moreover, in some of the above exam- 
ples a stress response crossover was reported either with 
system size and/or packing arrangement. This suggests 
that stress response investigations on granular packings 
can provide a route to understanding stress states in a 
range of systems where the discreteness of the constituent 
particles come into play. 



However, there currently remains a distinct lack of a 
systematic study on the role that system size plays in 
modifying the stress response in granular materials. We 
propose that such a study is needed with the ever-growing 
emphasis on the design of small devices where the par- 
ticulate nature of the material eventually dominate the 
system properties. Moreover, from this study we can es- 
timate relevant length scales over which continuum elas- 
ticity theory remains valid, and below which alternative 
descriptions may be needed. Thus, our goal here is to 
provide a qualitative guide map to the typical system 
parameters for which classical elasticity theory is likely 
to apply and those for which more advanced techniques 
may be needed to predict their properties. 



This paper is organized as follows. In the Section 
mi we summarize the Boussinesq equations of classi- 
cal, isotropic, linear elasticity theory for an infinite 
half-plane, and discuss other possibilities that predict 
anisotropic stress response profiles. In Section IIIII we 
provide an overview of the computational technique and 
the packing preparation protocol implemented in this 
work. Section HVl presents our main findings for the stress 
response behavior in frictional and frictionless particle 
packings. We end with conclusions and discussion. 



II. BOUSSINESQ EQUATIONS 

The equations of classical elasticity theory are derived 
from differential equations of the elliptic class that rely 
on the identification of well defined displacement fields 
of elements constituting the continuous elastic medium 
[l9| . The analytical solutions for the problem of a lo- 
calized point force applied vertically at the surface of an 
elastic, isotropic, infinite half-space are called the Boussi- 
nesq equations [l|, [HI. In Cartesian coordinates, where 
the origin of the coordinate system coincides with the ap- 
plication point of the vertical force, Fapp, the components 
of the displacement fields are given by. 
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where p = ^J ^ ^ . Both G, the modulus of rigid- 
ity, and the Poisson ratio, are material parameters 
characterizing the bulk properties of the elastic medium. 
Based on thermodynamic grounds [l9|, the Poisson ratio 
is constrained by, —1 < v < 0.5, in three dimensions. 
Through implementation of Hooke's law between stress 
and strain components, one then obtains expressions for 
components of the stress tensor. For the specific case 
when the force is applied in the 'downwards' or negative 
z-direction (see Fig. [Tj) the normal stress in the direction 
of the applied force, azz^ characterizes the response of 
the system. We focus primarily on this component of 
the stress although we also provide information on the 
normal stress components not in the direction of the ap- 
plied force, (J XX and (jyy. We give expressions for all three 
normal stress components 
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where r = x'^ + y^- 

The stress response scales linearly with the applied 
force. Deviations from this indicate non-linear behavior, 
which, in the case of a particulate medium, corresponds 
to non-affine, or irreversible, particle displacements in re- 
sponse to the applied load. We reiterate that that there 
are no fitting parameters in azz for the case of a semi- 
infinite half space. Therefore, we can quantify the influ- 
ence of boundaries, or confinement, through a compari- 
son of the Boussinesq result of Eq. [21 with our computed 
stress response for confined systems. Furthermore, we 
can also quantify deviations from isotropy from the stress 
profiles. We find that although the application of Boussi- 
nesq theory might at first seem to be inappropriate for 
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the systems investigated here, we find that for packings 
of sufficient linear size the classical theory captures the 
essential features of the stress properties and we classify 
the conditions for which isotropic elasticity theory are 
suitable. 

III. SIMULATION METHODS 

A. The Contact Forces and Equations of Motion 

We implement a granular dynamics (GD) variant on 
molecular dynamics simulations |20| that has been de- 
signed to simulate granular systems [2l|. Here we pro- 
vide a summary of the technique {llj. For this particu- 
lar study we focus on noncohesive, monodisperse sphere 
packings composed of elastic spheres of diameter d and 
mass m, that interact only on contact, through forces 
that act in directions that are normal n, and tangential 
t, to the contact plane, defined via 



t = 1 - n 

where Vij = — r^, is the separation between particles 
i and j, located at positions and vj respectively, and 
Tij = |r^j|. 1 represents the unit vector. 

The normal and tangential contact forces F^^t are 
based on a linear spring-dashpot model characterized by 
normal and tangential linear spring constants /Cn,t, and 
damping factors 7n,t, that account for elastic deformation 
of the contact point and inelasticity respectively, 

(4) 

Ft = -hsij - m^^-itv\.. 

Here, 5ij = {vij — d) is the surface compression of two 
particles undergoing a collision and v'^- the relative nor- 
mal/tangential velocity. The effective mass, meff = m/2, 
for the systems studied here. The quantity Sij represents 
the integrated displacement of the contact point while 
two particles remain in contact. Friction is implemented 
through a local Coulomb yield criterion: Ft < iiFn^ for a 
given friction coefficient fi. In this work, /i takes on the 
values zero and one. Thence, the total contact force is 
given by F = F^n + F^t. Interactions between the parti- 
cles and the walls are given by similar expressions where 
the wall is characterized by an effective mass of unity 
and the wall friction coefficient fiw takes on the same 
values as that of the particles. The results are presented 
in simulation units. Len gth scales are given in units of 
timescales in units of \/d/ where g is the acceleration 
due to gravity, forces in units of mg^ stresses in units of 
mg/d^^ and energies are given in units of mgd. Table [J 
shows the simulation parameters used in this study. We 
point out that the values of the spring constants are fixed 
throughout this study resulting in a particle scale Poisson 
ratio of zero. 



TABLE L Physical parameters used in this simulation study. 
kn and kt are normal and tangential spring constants in units 
of mg/d. 7n and are normal and tangential damping factors 
in units of \/ gj d. For values used here the particle Poisson 
ratio is zero. The particle friction coefficient is /i and the 
coefficient of inelasticity e. 



Number of Particles 


N 


1000 - 100000 


Normal Stiffness 


kn 


1 xlO^mg/d 


Tangential Stiffness 


kt 


1 xlO^mg/d 


Normal Damping Coefficient 


In 


50y^g/d 


Tangential Damping Coefficient 


In 


7n/2 


Normal Restitution Coefficient 


e 


0.88 


Particle Friction Coefficient 




0, 1.0 


Wall Friction Coefficient 




0, 1.0 



B. Coarse Graining Procedure 

Stresses are computed from the contact forces. To en- 
able an accurate determination of stresses from micro- 
scopic to macroscopic length scales, we follow the proto- 
col developed by Goldhirsch and colleagues [23]. We im- 
plement a coarse graining procedure into the microscopic 
stress expression [24|. The coarse graining function cho- 
sen is a 3D Gaussian, and the resulting expression for the 
stress components at some location r is, 

rr ^M--VF-- r--R t ds - e"^l^"^^+^^^^l^'/^' 

(5) 

where uj is the coarse graining length scale. Fij repre- 
sents the force experienced by particle i in contact with 
particle j, obtained from Eq. [H The numerical prefactor 
(V^cj)"^ is the normalization factor. 

To implement Eq. [5] into our calculations we adopted 
a spatial grid over which to evaluate the coarse grained 
stress. To achieve sufficient spatial resolution we fixed the 
grid size = O.ld and uj = l.Od. (We computed stresses for 
a range of uj and grid sizes for our 3D systems and found 
that O.bd < UJ < 5d provides a suitable plateau window 
[25I, at an acceptable computational cost.) 

C. Packing Generation Protocol 

In this work we investigate how confinement infiuences 
the nature of stress response in granular packings. Be- 
cause of the potentially large parameter space we fixed 
many features of the packings and focus our attention 
on a subset of variables. More specifically, we studied 
face centered cubic (FCC) ordered arrays of monodis- 
perse spheres confined within a walled box. To prepare 
the initial configuration, we placed non overlapping par- 
ticles at the lattice points of an FCC structure at the 
hard sphere packing fraction, 0fcc = 7r/\^Ts. The pack- 
ings are confined within open- top boxes, consisting of fiat 
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side walls and a bottom base. We defined a coordinate 
system as shown in Fig. [TJ The bottom plane of the box 
defines the x7/-plane, while gravity points in the —z direc- 
tion. The size of the simulation boxes were chosen to be 
commensurate with the crystal structure. For the FCC 
structures the lattice unit is defined as a = \f2d. Most of 
the results presented here are for boxes with square bases 
of side lengths, L, ranging from, 7a ^ lOd to 70a ^ lOOd. 
We also studied rectangular shaped boxes. The total 
number of particles A^, was chosen by fixing the height, 
iJ, of all packings to be, R ^ 10(i. 




FIG. 1: A computer simulation snapshot of a typical system 
studied in this work. The packing is composed of monodis- 
perse spheres arranged into a face centered cubic array, con- 
fined by 4 side walls, a solid base with an open top. The 
application of a localized perturbative force is represented by 
the arrow and the coordinate system is defined in the figure. 
X = and y — ^ define the center of box in the bottom of 
box. Gravity points in the —z direction. 

After the particles were positioned at the lattice sites, 
we switched on gravity and allowed the particles to settle 
until the energy per particle, EjN ^ 10~^^mgd^ reached 
numerical precision. This then defines our initial config- 
uration. We recorded the positions of the particles, the 
contact forces between interacting particles, and conse- 
quently the stress, ai, of this initial state. We then iden- 
tified the top middle region of the packing, locating the 
particles in this region. 

The second part of the simulation involves applying 
a downwards-pointing localized perturbation force, Fapp, 
to the top middle particle (s), and allowing the system 
to relax back to a mechanically stable state. We varied 
the magnitude of this applied force over several orders of 
magnitude, O.Olm^ < Fapp ^ lOOm^, to check the limits 
of the linear regime. Again, we tracked particle pos itions 
and contact forces, and the final stress state af [26|. The 
response function is calculated as the difference between 
the final and initial stress states, which we denote simply 



as, 

cr = CTf - (Ji. (6) 

To make a clearer connection with previous experimen- 
tal studies and visualizations [H, [l^ • we present our re- 
sults in a convenient format that can be readily compared 
with previous studies and theoretical predictions. In par- 
ticular, we make use of three dimensional contour stress 
maps and two dimensional stress profiles which view the 
stress response in a given plane of the box or angularly- 
averaged slice across this plane, respectively. Given the 
box geometry that we employ, it is convenient to visualize 
(Jzz as a stress map in the horizontal plane. Using these 
plots it is relatively straightforward to distinguish be- 
tween isotropic and anisotropic response profiles. These 
typically show up more clearly in stress profiles as single- 
and double-peak response profiles respectively. In all our 
stress maps lighter shading indicates larger stresses. We 
also present results on particle displacement fields taken 
in the plane. 

IV. RESULTS 

A. Frictional Packings 

Our first results show that ordered arrays of frictional 
granular packings display isotropic elastic-like stress re- 
sponse behavior as seen in Fig. [2l These stress maps 
convey the "single-peak" nature of the stress response 
measured at the bottom of the packing for different sys- 
tem sizes for one value of the applied force, Fapp = Img. 
Thus, we immediately note that the Boussinesq equa- 
tion appears to describe confined systems. However, we 
find that system size and force magnitude influence this 
picture. We quantify the regime over which the classical 
expression remains suitably valid in two ways: Firstly, we 
compare our stress profiles to the Boussinesq expression 
and extract the full width at half maximum. Secondly, 
we identify the linear regime as a function of the applied 
force. 




FIG. 2: Frictional systems /i = 1. Stress response maps of azz 
at the bottom boundary due to an applied force Fapp = Img, 
for different system sizes, lOd < L < lOOd indicated by scale 
in panels. 

To cast this data into a more familiar form we circu- 
larly average the stress map data of Fig. [2] to construct 
stress profiles and directly compare with the predictions 
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of Eq. [21 from which we also obtain the fuh width of 
the response peaks. In Fig. [31(a) we show our averaged 
response profiles for systems of different height, clearly 
indicating a single-peak response consistent with the clas- 
sical theory shown by the thick solid line. We note that 




FIG. 3: (a) Averaged stress response profiles, a, for one value 
of the applied force Fapp = Img for systems of fixed height 
H ^ lOd and varying side lengths L = lOd (long dash line), 
20d (thin solid line), 40d (dashed), lOOd (dotted), and B (solid 
line), the Boussinesq result, (b) A full width at half height of 
the stress profiles for the simulation data (symbols) relative 
to the theoretical result Ab- Perfect agreement indicated by 
dashed line. 

for the smallest system {L = lOd) the stress profiles do 
not accurately conform to the Boussinesq form. To bet- 
ter quantify deviations from Eq. [2] we extract the full 
width at half height, A, and compare with the theoret- 
ical result A^ = 1.113i^, shown in Fig. [3^b). We find 
that for L > 40(i, the single-peak response profiles are 
suitably described by Eq. [21 with better agreement for 
larger system size as expected. 

Thus, we have identified a minimum system size for 
which the classical expression describes the data accu- 
rately, L > AOd. To further test our results with Eq. [21 
we again fit the data and now allow the applied force pa- 
rameter in Eq. [2]to vary as a fitting parameter, denoted 
here as Fg, which we then compared to the actual force 
applied during the simulation, Fapp. Our results shown 
in Fig. [4] clearly indicate that confined systems of box 
sizes greater than L > AOd conform to the semi-infinite 
system size result. 

1 .5 I — I — I 1 1 — 

1 

0.5 - * 

I ^ ^ ^ ^ 

10 20 40 100 

L/d 

FIG. 4: Comparison between the actual applied force Fapp — 
Img versus Fb, the fitted force using the Boussinesq equation, 
Eq. [21 to fit the stress profile data at the bottom of a frictional 
granular packing for boxes of different sizes, with side lengths 
lOd < L < lOOd. Perfect agreement indicated by dashed line. 

The magnitude of the applied force also strongly influ- 



ences the response function. For small applied forces, the 
stress response remains isotropic consistent with the clas- 
sical result as seen in the previous figures. Beyond some 
force threshold the stress response becomes anisotropic 
due to the FCC arrangement of the grains and Eq. [2] is 
no longer valid as indicated by the stress maps in Fig. [5l 
To quantify this effect and identify the linear regime as 
a function of applied force, we performed a similar anal- 
ysis as described above. We varied the applied force over 
several orders of magnitude, O.Olm^ < Fapp < lOOm^, 
applied to a system size of side length L = 40(i. Figure 
[6] shows that the fitted force increases linearly with ap- 
plied force with unit slope when the magnitude of the 
perturbation force remains below Fapp — 60mg. In ongo- 
ing work to be reported elsewhere [28] we find that this 
threshold also depends on the friction coefficient. Thus, 




FIG. 5: Stress response map, cr^^, for the bottom of the box of 
a frictional system of size L = 20d, Left panel: Fapp = lOm^. 
Right panel: Fapp = lOOm^'. 
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FIG. 6: Comparison between the actual applied force Fapp 
versus Fb, the fitted force using the Boussinesq equation to 
fit the stress profile data at the bottom of a frictional granular 
packing of side length L = 20d. 

in summary, we find that for ordered frictional pack- 
ings the stress state is suitably described by the semi- 
infinite result even for relatively small systems over a 
wide range in applied forcing and we have quantified the 
linear regime over which this agreement exists. 

In an effort to investigate the underlying mechanisms 
responsible for the changing character of the stress re- 
sponse in our different systems we have computed the 
normal stress response at different distances from the 
source of the force perturbation which are shown in 
Fig. [3 The left hand panels in Fig. [7] show data for 
our smallest system, L = lOd, which deviates strongly 
from the classical predictions despite appearing to retain 
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a single broad peak response profile. While the right 
hand columns are for a larger system, L = 40(i, that is 
consistent with the Boussinesq expression. As expected 
at the top of packing (top panels) the stress is localized 
at the point of perturbation. Deeper into the packings 
the stress response broadens. For the smallest system 
the stress quickly spreads across the entire region of the 
packing while for larger systems the stress peak remains 
relatively localized within the central region of the pack- 
ing. 




0.02 
0.01 



FIG. 7: Stress maps, azz, in response to a localized force of 
magnitude Fapp = Img applied at the top of the packing. The 
panels show data at different distances from the perturbation 
source. Rows show data for layer 1 (top), 3, and 8 (near 
bottom), for L = lOd (left column panels) and L — 4Qd (right 
column panels). 

The reason behind this changing stress response is a 
strong reflective component of stress from the side walls 
as the stress is transmitted away from the source of the 
perturbation. To qualify this effect, in Fig. [8] we plot the 
averaged normal stresses, \{cFxx ^yy) = '''^xx \ exerted 
at the side walls for our different systems, and data for 
the classical expression using a value of the Poisson ratio 
of = 0.4. Indeed, we find that for the smaller systems 
the stress at the side walls is much larger than that ex- 
pected from the classical behavior of side wall stress from 
Boussinesq equation and the corresponding larger system 
size. To generate the data presented in Fig. [8] it is neces- 
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FIG. 8: Normal stress, Oxx^ at the side wall for frictional 
packings of height ^ lOc? and side length, L = lOd, 20d, 40d 
and the classical expression for a system of size L — 40(i 
(bottom right panel), using v — 0.4. 

sary to introduce the material Poisson ratio v. Our final 
choice for z/ was made on the best fit between the stress 



profiles obtained via the displacement fields, Eq. [H and 
the displacement fields constructed from the simulation 
data. We matched our simulation displacement fields for 
our largest system to the predictions of Eq. [1] using u 
as a fitting parameter. Our confined packings therefore 
behave as a continuum material with an effective Pois- 
son ratio, V ~ 0.4. This procedure actually provides a 
convenient method to determine material properties of 
composite materials such as granular packings [29|. The 
benefit of analyzing the displacement fields is that we 
can view the direct microscopic response of the system 
at the particle scale due to the imposed localized force 
perturbation. 

The displacement fields of the middle vertical slice for 
frictional packings of different size are plotted in Fig. [9l 
We also include a comparison to Eq. [T] for a system of 
size L = 100, using u = 0.4. These changes in the dis- 
placement fields of the perturbed system result in a larger 
reflection of stress at the side walls hence causing devi- 
ations in the overall stress response of smaller systems 
compared with the semi-infinite size result. However, it 
is worth pointing out that these differences occur only for 
relatively small packings of size L < 20d. Above this size, 
these frictional packings appear to be suitably described 
by the Boussinesq formalism. 
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FIG. 9: Displacement field vectors in the vertical plane for a 
central slice inside packings of side length L = 10, 20, lOOd, 
and the result of Eq. [T] on a system of size L = lOOd, using 
ly = 0.4. 



B. Frictionless Packings 

In contrast to frictional packings, frictionless ordered 
arrays display primarily strongly anisotropic stress be- 
havior in response to localized force perturbations. In 
Fig. [To] the stress component cizz for frictionless arrays of 
different sizes are compared. Apart from the smallest sys- 
tem size, the response profiles are ringed or multi-peaked 
with a minimum in the stress response characteristic of 
an anisotropic response function in clear contrast to the 
frictional packings of Fig. [2l Surprisingly, however, the 
smallest system of size L = lOd exhibits a stress response 
that initially appears isotropic in nature - a maximum in 
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the stress in the middle of the packing - but with some 
unusual features not seen in the frictional case. In fact 
this smaller system does not conform to the Boussinesq 
stress profile even though an isotropic-like single peak 
response is observed. 




FIG. 10: Frictionless systems /i = 0. Stress response maps of 
Gzz at the bottom boundary due to an applied force 
Img, for different systems sizes, lOd < L < lOOd. 

The influence of confinement plays a significant role 
in modifying the stress properties for frictionless pack- 
ings even within the same packing. To investigate this 
effect further, stress response maps at different distances 
from the source of the force perturbation are shown in 
Fig. [11] for two system sizes, L = lOd, 40(i. In the vicin- 
ity of the localized force perturbation, the stress response 
also takes on localized features where the stress is con- 
centrated in a region directly below the point of appli- 
cation. At intermediate distances between the top and 
bottom of the packing both systems now display strongly 
anisotropic stress profiles with a stress minimum focused 
in the region beneath the point of force application. How- 
ever, as noted above, the response measured at the bot- 
tom of the packing dramatically changes character be- 
tween the two system sizes. These features suggest that 
the stress state of small frictionless packings can not only 
vary substantially within the packing, but exhibit highly 
unusual stress properties compared to both larger fric- 
tionless systems and frictional packings. 




FIG. 11: Stress maps, azz, in response to a localized force of 
magnitude Fapp = Img applied at the top of the packing. The 
panels show data at different distances from the perturbation 
source. Rows show data for layer 1 (top), 3, and 8 (near 
bottom), for L = lOd (left column panels) and L = 40d (right 
column panels). 

Again we find that these unusual stress properties are 



a direct result of the manner in which stresses are trans- 
mitted and interact with the side wall boundaries due 
to the strongly correlated displacement fields of the par- 
ticles during the response process. Visualization of the 
displacement fields for these frictionless packings are par- 
ticularly illuminating as shown in Fig. [121 These dis- 
placement vectors are taken from a middle slice inside 
the different packings. For the smallest packing (top left 
panel in Fig. [12]), the particle displacements in response 
to the applied perturbation show a directional structure 
that is responsible for the increased stress reflection at 
the side walls. Interestingly, this reflection shields the 
bottom of the packing from the growing stress minimum 
apparent at intermediate depths. Whereas, for larger 
systems the stress response is seen to be dominated by 
particle displacements along directed rays through the 
system responsible for the anisotropic profiles observed 
at the bottom of the packings. 
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FIG. 12: Displacement field vectors in the vertical plane for a 
central slice inside packings of side length L = 10, 20, 40, lOOd. 

The results presented above suggest that there exists 
some crossover response for confined, frictionless pack- 
ings that depends on the size of the system. We have in- 
vestigated the changing character of the stress response 
and its dependence on box size. This crossover behavior 
is illustrated in stress response maps of Fig. [131 where 
the box side lengths have been incrementally increased 
from L = 12d — 18d. It is apparent that changing box 




FIG. 13: Crossover behavior of the stress response, cr^z, with 
increasing box size for systems of side length, 12d < L < 18d. 

size results in large-scale stress fluctuations that span 
the system as the stress response crosses over from a 
more isotropic-like response (smaller systems) to strongly 
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anisotropic behavior (larger system) . There is a region of 
box sizes for these frictionless packings where the stress 
state of the system exhibits pecuhar behavior that does 
not conform to either a pseudo-isotropic description nor 
can it be considered anisotropic in character. These fea- 
tures can have imphcations in the design of mechanical 
systems at small scales where confinement can have a 
large effect on the stress state of the system. 

We have also investigated the role of box geometry on 
the nature of the stress response. As shown earlier we find 
that for the smaller frictionless systems confined within a 
square-base box with = Ly = lOd, the stress response 
exhibits a broad peak in the central region of the base. 
We find that highly unusual and complex features emerge 
if we now increase the size of only one side of the box 
into a rectangular base keeping Lx = lOd fixed. Stress 
maps for this scenario are shown in Fig. [141 where we 
have increased the side ratio from 1 to 10. These result 
emphasize the importance of box size and geometry on 
the nature of stress response in shallow packings. 
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FIG. 14: Evolution of the stress response, cr^^, in frictionless 
packings as one side of the confining box is increased from a 



square base of side lengths Lx 
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V. CONCLUSIONS 

Using computer simulations of a model granular mate- 
rial, we have studied the nature of the stress response in 
shallow and confined granular crystalline arrays. Pack- 
ings composed of highly frictional particles generally ex- 
hibit a stress response function measured at the bottom 
of the packing that resembles that of an isotropic, elastic 
material that is well described by the Boussinesq equa- 
tions of classical elasticity. Thus, the semi-infinite half- 
space result remains suitably valid to describe the stress 
state of the system provided the size of the base of the 
confining walls is larger than approximately 40 particle 
diameters. In this regime the stress response appears to 
behave as a linear elastic material. For smaller frictional 
systems, the stress response exhibits confinement effects 
that are not accurately captured by the classical result. 
Packing arrangement and the infiuence of the boundaries 
have significant effects on the stress state of the system. 
Furthermore, due to the FCC crystalline arrangements 



studied here, anisotropic stress response behavior is ob- 
served at large applied forces whereby the infiuence of 
structure starts to dominate the stress properties. 

Frictionless packings exhibit stress response proper- 
ties that are generally anisotropic in nature due to the 
FCC arrays. However, smaller packings can be tuned to 
observe a range of behaviors depending on the location 
within the packing for square boxes, or on the ratio of the 
side lengths for rectangular systems. Complex response 
patterns emerge for frictionless rectangular arrays. We 
also point out that for the anisotropic profiles - either 
with zero friction or for large forcing - the resulting stress 
response indicates a decrease in the stress relative to the 
initial state prior to perturbation, i.e. in the measures 
defined here the stress becomes negative. This suggests 
that such force perturbations actually weaken the mate- 
rial relative to the unperturbed state. 

Our results further suggest that there exists a transi- 
tion in the underlying character of the response as the 
friction coefficient is varied. To highlight this possibil- 
ity, we show in Fig. [151 the averaged stress profiles for 
a square-base, L = lOOd, for different particle friction 
coefficients in response to an applied force, Fapp = Img. 
Indeed we find that increasing the friction coefficient sup- 
presses the stress dip in the central region of the packing 
until a single-peak response is recovered for moderate 
to large friction coefficients. We expect that the precise 
location of this transition will also depend on the magni- 
tude of the applied force psj . 




FIG. 15: cr, averaged stress profiles for different friction coef- 
ficients, /i, for one system size, L = lOOd, in response to an 
applied force Fapp = Img. Each stress profile is scaled by its 
maximum value, am 



Our studies indicate that for confined granular arrays, 
particle properties, magnitude of forcing, and geometry 
of the confining container can all be used to design mate- 
rials with specific stress states. Through this accessible 
parameter space, confined granular packings can be made 
to mimic the properties of linear, isotropic, elastic mate- 
rials, where the stress is concentrated beneath the point 
of perturbation - single-peak response. Or they can be 
designed to exhibit strongly anisotropic properties where 
stresses are transmitted along preferred directions. More- 
over, for particular box geometries, the stress response 
can be transmitted approximately uniformly through the 
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system. Thus, it might be necessary to take into account 
the possible stress states indicated here during the design 
of microscopic granular devices and components. 



VI. APPENDIX 

Here we provide a short derivation of the normal stress 
response component a^z (see Eq. [2|), within the frame- 
work of linear elasticity, indicating how the material pa- 
rameters, G and cancel from the expression given in 
Eq. O This derivation follows directly from Landau and 
Lifshitz [l9|. The basic ingredient is the implementa- 
tion of a constitutive relation between stress and strain 
that corresponds to Hooke's law: the normal stress a^z is 
linearly proportional to the strain components u^x^ Uyy, 



and Uzz^ 



O'zz = 



2G 
-2v 



(7) 



where the strain components are obtained from deriva- 
tives of the components of the displacement fields leading 
to 



"^xx — 



47rG 



47rG 



3a; z 



p5 



l-2v 



l-2v 



{p^+pzY 



{l-2v)y\2+^^ 



{p^+pzY 



(8) 



1 

47rG 



2vz 

n3 



and we have set the value of the applied, perturbing force 
to unity for notational convenience. 



Combining the above equations we find. 



2(1-2^) (l-2^)(x2+^2)(2 + ^) 2vz ?>z\ 



pz 



{p^^pzf 



explicitly emphasizing the cancellation of G from the 
equation. 

The resulting expression can be written as 



2^1 -2Z.) Kl-2^)(l-f)(2 + f)^ 



- pz 



and with further algebra. 



pz 



CTz 



1 

2^ 



vz 32: 1 , uz uz\ 

p'^ p^ p"^ -\- rz p 



The terms containing u cancel out and one obtains Eq. [2] 
for an arbitrary applied force, Fapp, directed in the —z 
direction. 
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